probabilistic programming

modelling workflow for risk and safety

Dom Di Francesco

prepared for the Health & Safety Executive, Science Division

13 Feb 2025

overview

introduction to PPLs
modelling workflow:
- prior checks
- sampling diagnostics
- model evaluation
- experimental design
…if we have time:
- measurement error
- mixture models
- hierarchical models

from `BUGS`/`JAGS` to `Stan`

Gibbs sampling:
- draw samples from one parameter at a time
- helpful when Metropolis-Hastings is rejecting proposals
Hamiltonian Monte Carlo:
- uses gradient of log-posterior and Hamiltonian dynamics analogy
- samples from all continuous parameters jointly and efficiently
- extended to No-U-Turn Sampler (NUTS) in Stan, which tunes hyperparameters during warmup!
link to a fun demo

from `BUGS`/`JAGS` to `Stan`

Hoffman, M.D., Gelman, A., ‘The No-U-Turn Sampler’, JMLR, 2014

from `BUGS`/`JAGS` to `Stan`

StanCon anyone?

PPLs

‘democratising scalable UQ’ …once installed 🤔

Stan (via CmdStanR)
Stan (via CmdStanPy)
Turing.jl

Code

library(cmdstanr)

cmdstanr::cmdstan_version()

[1] "2.34.1"

Code

import cmdstanpy

cmdstanpy.cmdstan_version()

(2, 34)

Code

using Turing, Pkg

Pkg.status("Turing")

Status `~/.julia/environments/v1.11/Project.toml`
⌃ [fce5fe82] Turing v0.35.4
Info Packages marked with ⌃ have new versions available and may be upgradable.

corrosion rates

looking for evidence of active corrosion growth

corrosion rates

Code

library(tidyverse)

corrosion_data <- read_csv("../../data/corrosion.csv")
head(corrosion_data, 3)

# A tibble: 3 × 5
  anomaly_id soil_type inspection measured_depth_mm     T
       <dbl> <chr>     <chr>                  <dbl> <dbl>
1          1 A         i_0                   0.0672     0
2          2 A         i_0                   0.104      0
3          3 A         i_0                   0.101      0

Code

import polars as pl

corrosion_data = pl.read_csv("../../data/corrosion.csv")
corrosion_data.head(3)

shape: (3, 5)

anomaly_id	soil_type	inspection	measured_depth_mm	T
i64	str	str	f64	i64
1	"A"	"i_0"	0.067209	0
2	"A"	"i_0"	0.104202	0
3	"A"	"i_0"	0.100588	0

Code

using CSV, DataFrames

corrosion_data = CSV.read("../../data/corrosion.csv", DataFrame)
first(corrosion_data, 3)

3×5 DataFrame
 Row │ anomaly_id  soil_type  inspection  measured_depth_mm  T
     │ Int64       String1    String3     Float64            Int64
─────┼─────────────────────────────────────────────────────────────
   1 │          1  A          i_0                 0.0672091      0
   2 │          2  A          i_0                 0.104202       0
   3 │          3  A          i_0                 0.100588       0

a corrosion growth rate model

\[\begin{aligned} & \Delta C = \frac{C_{j} - C_{i}}{\Delta t_{i \rightarrow j}} \\ \\ & \Delta C \sim N(\mu, \sigma^2) \end{aligned}\]

CmdStanR
CmdStanPy
Turing.jl

Code

cgr_model <- cmdstan_model(stan_file = "corrosion_growth.stan")

cgr_model$format()

data {
  int<lower=0> n_anomalies; // number of anomalies
  int<lower=0> n_inspections; // number of inspections
  vector[n_anomalies * n_inspections] cgr; // growth rate observations (2 per anomaly)
}
parameters {
  real<lower=0> mu; // mean growth rate
  real<lower=0> sigma; // standard deviation
}
model {
  // model
  for (i in 1 : (n_anomalies * n_inspections)) {
    cgr[i] ~ normal(mu, sigma);
  }
  /*
  //alternative (vectorised) implementation:
  delta_C ~ normal(mu, sigma);
  
  //some suggested priors
  mu ~ normal(1/4, 3);
  sigma ~ exponential(1);
  */
}
generated quantities {
  vector[n_anomalies * n_inspections] cgr_pred;
  //vector[n_anomalies * n_inspections] log_lik;
  
  for (i in 1 : (n_anomalies * n_inspections)) {
    cgr_pred[i] = normal_rng(mu, sigma);
    //log_lik[i] = normal_lpdf(delta_C[i] | mu, sigma);
  }
}

Code

cgr_model = cmdstanpy.CmdStanModel(stan_file="corrosion_growth.stan")

INFO:cmdstanpy:found newer exe file, not recompiling

Code

stan_code = cgr_model.code()

from pygments import highlight
from pygments.lexers import StanLexer
from pygments.formatters import NullFormatter

formatted_stan_code = highlight(stan_code, StanLexer(), NullFormatter())

print(formatted_stan_code)

data {
  int <lower = 0> n_anomalies;  // number of anomalies
  int <lower = 0> n_inspections;  // number of inspections
  vector[n_anomalies * n_inspections] cgr;  // growth rate observations (2 per anomaly)
}

parameters {
  real<lower = 0> mu;  // mean growth rate
  real<lower = 0> sigma;  // standard deviation
}

model {  
  // model
  for(i in 1:(n_anomalies * n_inspections)){
    cgr[i] ~ normal(mu, sigma);
  }
  /*
  //alternative (vectorised) implementation:
  delta_C ~ normal(mu, sigma);

  //some suggested priors
  mu ~ normal(1/4, 3);
  sigma ~ exponential(1);
  */
  
}

generated quantities {
  vector[n_anomalies * n_inspections] cgr_pred;
  //vector[n_anomalies * n_inspections] log_lik;

  for (i in 1:(n_anomalies * n_inspections)) {
    cgr_pred[i] = normal_rng(mu, sigma);
    //log_lik[i] = normal_lpdf(delta_C[i] | mu, sigma);
  }
}

Code

@model function corrosion_growth(cgr)
    # priors
    μ ~ Normal(0, 2) |> d -> truncated(d, lower = 0)
    σ ~ Exponential(1)
    
    # model
    for i in eachindex(cgr)
        cgr[i] ~ Normal(μ, σ) |> d -> truncated(d, lower = 0)
    end

    # Turing automatically keeps track of log-likelihoods 🏆 
    
end

running the model

CmdStanR
CmdStanPy
Turing.jl

CmdStanR needs it’s input data as a list

Code

prepare_data <- function(df = corrosion_data) {
  df |>
    arrange(anomaly_id, T) |>
    group_by(anomaly_id) |>
    mutate(
      next_depth = lead(measured_depth_mm),
      time_diff = lead(T) - T
    ) |>
    filter(!is.na(next_depth)) |>
    mutate(
      delta_C = (next_depth - measured_depth_mm) / time_diff
    ) |>
    select(anomaly_id, delta_C, soil_type) |>
    ungroup()
}

model_data <- list(
  n_anomalies = prepare_data()$anomaly_id |> unique() |> length(),
  n_inspections = 2,
  cgr = prepare_data()$delta_C
)

cgr_post <- cgr_model$sample(data = model_data)

Running MCMC with 4 sequential chains...

Chain 1 Iteration:    1 / 2000 [  0%]  (Warmup) 
Chain 1 Iteration:  100 / 2000 [  5%]  (Warmup) 
Chain 1 Iteration:  200 / 2000 [ 10%]  (Warmup) 
Chain 1 Iteration:  300 / 2000 [ 15%]  (Warmup) 
Chain 1 Iteration:  400 / 2000 [ 20%]  (Warmup) 
Chain 1 Iteration:  500 / 2000 [ 25%]  (Warmup) 
Chain 1 Iteration:  600 / 2000 [ 30%]  (Warmup) 
Chain 1 Iteration:  700 / 2000 [ 35%]  (Warmup) 
Chain 1 Iteration:  800 / 2000 [ 40%]  (Warmup) 
Chain 1 Iteration:  900 / 2000 [ 45%]  (Warmup) 
Chain 1 Iteration: 1000 / 2000 [ 50%]  (Warmup) 
Chain 1 Iteration: 1001 / 2000 [ 50%]  (Sampling) 
Chain 1 Iteration: 1100 / 2000 [ 55%]  (Sampling) 
Chain 1 Iteration: 1200 / 2000 [ 60%]  (Sampling) 
Chain 1 Iteration: 1300 / 2000 [ 65%]  (Sampling) 
Chain 1 Iteration: 1400 / 2000 [ 70%]  (Sampling) 
Chain 1 Iteration: 1500 / 2000 [ 75%]  (Sampling) 
Chain 1 Iteration: 1600 / 2000 [ 80%]  (Sampling) 
Chain 1 Iteration: 1700 / 2000 [ 85%]  (Sampling) 
Chain 1 Iteration: 1800 / 2000 [ 90%]  (Sampling) 
Chain 1 Iteration: 1900 / 2000 [ 95%]  (Sampling) 
Chain 1 Iteration: 2000 / 2000 [100%]  (Sampling) 
Chain 1 finished in 0.1 seconds.
Chain 2 Iteration:    1 / 2000 [  0%]  (Warmup) 
Chain 2 Iteration:  100 / 2000 [  5%]  (Warmup) 
Chain 2 Iteration:  200 / 2000 [ 10%]  (Warmup) 
Chain 2 Iteration:  300 / 2000 [ 15%]  (Warmup) 
Chain 2 Iteration:  400 / 2000 [ 20%]  (Warmup) 
Chain 2 Iteration:  500 / 2000 [ 25%]  (Warmup) 
Chain 2 Iteration:  600 / 2000 [ 30%]  (Warmup) 
Chain 2 Iteration:  700 / 2000 [ 35%]  (Warmup) 
Chain 2 Iteration:  800 / 2000 [ 40%]  (Warmup) 
Chain 2 Iteration:  900 / 2000 [ 45%]  (Warmup) 
Chain 2 Iteration: 1000 / 2000 [ 50%]  (Warmup) 
Chain 2 Iteration: 1001 / 2000 [ 50%]  (Sampling) 
Chain 2 Iteration: 1100 / 2000 [ 55%]  (Sampling) 
Chain 2 Iteration: 1200 / 2000 [ 60%]  (Sampling) 
Chain 2 Iteration: 1300 / 2000 [ 65%]  (Sampling) 
Chain 2 Iteration: 1400 / 2000 [ 70%]  (Sampling) 
Chain 2 Iteration: 1500 / 2000 [ 75%]  (Sampling) 
Chain 2 Iteration: 1600 / 2000 [ 80%]  (Sampling) 
Chain 2 Iteration: 1700 / 2000 [ 85%]  (Sampling) 
Chain 2 Iteration: 1800 / 2000 [ 90%]  (Sampling) 
Chain 2 Iteration: 1900 / 2000 [ 95%]  (Sampling) 
Chain 2 Iteration: 2000 / 2000 [100%]  (Sampling) 
Chain 2 finished in 0.1 seconds.
Chain 3 Iteration:    1 / 2000 [  0%]  (Warmup) 
Chain 3 Iteration:  100 / 2000 [  5%]  (Warmup) 
Chain 3 Iteration:  200 / 2000 [ 10%]  (Warmup) 
Chain 3 Iteration:  300 / 2000 [ 15%]  (Warmup) 
Chain 3 Iteration:  400 / 2000 [ 20%]  (Warmup) 
Chain 3 Iteration:  500 / 2000 [ 25%]  (Warmup) 
Chain 3 Iteration:  600 / 2000 [ 30%]  (Warmup) 
Chain 3 Iteration:  700 / 2000 [ 35%]  (Warmup) 
Chain 3 Iteration:  800 / 2000 [ 40%]  (Warmup) 
Chain 3 Iteration:  900 / 2000 [ 45%]  (Warmup) 
Chain 3 Iteration: 1000 / 2000 [ 50%]  (Warmup) 
Chain 3 Iteration: 1001 / 2000 [ 50%]  (Sampling) 
Chain 3 Iteration: 1100 / 2000 [ 55%]  (Sampling) 
Chain 3 Iteration: 1200 / 2000 [ 60%]  (Sampling) 
Chain 3 Iteration: 1300 / 2000 [ 65%]  (Sampling) 
Chain 3 Iteration: 1400 / 2000 [ 70%]  (Sampling) 
Chain 3 Iteration: 1500 / 2000 [ 75%]  (Sampling) 
Chain 3 Iteration: 1600 / 2000 [ 80%]  (Sampling) 
Chain 3 Iteration: 1700 / 2000 [ 85%]  (Sampling) 
Chain 3 Iteration: 1800 / 2000 [ 90%]  (Sampling) 
Chain 3 Iteration: 1900 / 2000 [ 95%]  (Sampling) 
Chain 3 Iteration: 2000 / 2000 [100%]  (Sampling) 
Chain 3 finished in 0.1 seconds.
Chain 4 Iteration:    1 / 2000 [  0%]  (Warmup) 
Chain 4 Iteration:  100 / 2000 [  5%]  (Warmup) 
Chain 4 Iteration:  200 / 2000 [ 10%]  (Warmup) 
Chain 4 Iteration:  300 / 2000 [ 15%]  (Warmup) 
Chain 4 Iteration:  400 / 2000 [ 20%]  (Warmup) 
Chain 4 Iteration:  500 / 2000 [ 25%]  (Warmup) 
Chain 4 Iteration:  600 / 2000 [ 30%]  (Warmup) 
Chain 4 Iteration:  700 / 2000 [ 35%]  (Warmup) 
Chain 4 Iteration:  800 / 2000 [ 40%]  (Warmup) 
Chain 4 Iteration:  900 / 2000 [ 45%]  (Warmup) 
Chain 4 Iteration: 1000 / 2000 [ 50%]  (Warmup) 
Chain 4 Iteration: 1001 / 2000 [ 50%]  (Sampling) 
Chain 4 Iteration: 1100 / 2000 [ 55%]  (Sampling) 
Chain 4 Iteration: 1200 / 2000 [ 60%]  (Sampling) 
Chain 4 Iteration: 1300 / 2000 [ 65%]  (Sampling) 
Chain 4 Iteration: 1400 / 2000 [ 70%]  (Sampling) 
Chain 4 Iteration: 1500 / 2000 [ 75%]  (Sampling) 
Chain 4 Iteration: 1600 / 2000 [ 80%]  (Sampling) 
Chain 4 Iteration: 1700 / 2000 [ 85%]  (Sampling) 
Chain 4 Iteration: 1800 / 2000 [ 90%]  (Sampling) 
Chain 4 Iteration: 1900 / 2000 [ 95%]  (Sampling) 
Chain 4 Iteration: 2000 / 2000 [100%]  (Sampling) 
Chain 4 finished in 0.1 seconds.

All 4 chains finished successfully.
Mean chain execution time: 0.1 seconds.
Total execution time: 0.6 seconds.

CmdStanPy needs it’s input data as a dictionary

Code

def prepare_data(df = corrosion_data):
    return (
        df.sort(['anomaly_id', 'T'])
        .group_by('anomaly_id')
        .agg([
            pl.col('measured_depth_mm').shift(-1).alias('next_depth'),
            pl.col('T').shift(-1).alias('next_time'),
            pl.col('measured_depth_mm'),
            pl.col('T')
        ])
        .filter(pl.col('next_depth').is_not_null())
        .with_columns([
            ((pl.col('next_depth') - pl.col('measured_depth_mm')) / 
             (pl.col('next_time') - pl.col('T'))).alias('delta_C')
        ])
        .select(['anomaly_id', 'delta_C'])
        .explode('delta_C')  # Add this line to unnest the lists
        .filter(pl.col('delta_C').is_not_null())  # Optional: remove null values if any
    )

model_data = {
        'n_anomalies': prepare_data().select('anomaly_id').unique().height,
        'n_inspections': 2,
        'cgr': prepare_data().select('delta_C').to_series().to_numpy()
    }

cgr_post = cgr_model.sample(data = model_data)

                                                                                                                                                                                                                                                                                                                                

INFO:cmdstanpy:CmdStan start processing

chain 1 |[33m          [0m| 00:00 Status

chain 2 |[33m          [0m| 00:00 Status[A


chain 3 |[33m          [0m| 00:00 Status[A[A



chain 4 |[33m          [0m| 00:00 Status[A[A[A

chain 2 |[34m######### [0m| 00:00 Iteration: 1700 / 2000 [ 85%]  (Sampling)[A
chain 1 |[34m#########5[0m| 00:00 Iteration: 1800 / 2000 [ 90%]  (Sampling)


chain 3 |[34m########1 [0m| 00:00 Iteration: 1500 / 2000 [ 75%]  (Sampling)[A[A



chain 4 |[34m#########5[0m| 00:00 Iteration: 1800 / 2000 [ 90%]  (Sampling)[A[A[A
chain 1 |[34m##########[0m| 00:00 Sampling completed                       

chain 2 |[34m##########[0m| 00:00 Sampling completed                       

chain 3 |[34m##########[0m| 00:00 Sampling completed                       

chain 4 |[34m##########[0m| 00:00 Sampling completed                       
INFO:cmdstanpy:CmdStan done processing.

A Turing model needs it’s input data as arguments to the model function

Code

function prepare_data(df::DataFrame = corrosion_data)
    sorted_df = sort(df, [:anomaly_id, :T]); result = DataFrame()
    
    for group in groupby(sorted_df, :anomaly_id)
        if nrow(group) > 1
            for i in 1:(nrow(group)-1)
                Δc = (group[i+1, :measured_depth_mm] - group[i, :measured_depth_mm]) / (group[i+1, :T] - group[i, :T])
                push!(result, (
                    anomaly_id = group[i, :anomaly_id], Δc = max(0, Δc)
                ))
            end
        end
    end
    
    return result
end

cgr_post = prepare_data().Δc |> 
  data -> corrosion_growth(data) |>
  model -> sample(model, NUTS(), 1_000)

taking a look

R
Python
Julia

Code

?cmdstanr::draws()

Code

cgr_post$draws(format = "df")

# A draws_df: 1000 iterations, 4 chains, and 53 variables
   lp__    mu sigma cgr_pred[1] cgr_pred[2] cgr_pred[3] cgr_pred[4] cgr_pred[5]
1   130 0.083 0.046       0.013       0.085       0.070       0.097       0.114
2   129 0.087 0.046       0.083       0.100       0.039       0.054       0.051
3   130 0.079 0.041       0.064       0.041      -0.021       0.124       0.130
4   130 0.081 0.043       0.029       0.089       0.093       0.135       0.095
5   131 0.083 0.041       0.074       0.062       0.060       0.069       0.044
6   130 0.082 0.038       0.108      -0.016       0.111       0.064       0.094
7   130 0.079 0.043       0.087       0.049      -0.019       0.057       0.117
8   128 0.095 0.039       0.086       0.176       0.110       0.089       0.102
9   128 0.069 0.044       0.047       0.056      -0.011       0.127       0.088
10  129 0.071 0.040       0.090       0.059       0.027       0.109       0.035
# ... with 3990 more draws, and 45 more variables
# ... hidden reserved variables {'.chain', '.iteration', '.draw'}

Code

DomDF::tidy_mcmc_draws(cgr_post, params = c("mu", "sigma"))

# A tibble: 8,000 × 4
   Parameter Chain Iteration  value
   <chr>     <int>     <int>  <dbl>
 1 mu            1         1 0.0830
 2 mu            1         2 0.0875
 3 mu            1         3 0.0792
 4 mu            1         4 0.0807
 5 mu            1         5 0.0827
 6 mu            1         6 0.0824
 7 mu            1         7 0.0793
 8 mu            1         8 0.0948
 9 mu            1         9 0.0690
10 mu            1        10 0.0712
# ℹ 7,990 more rows

Code

cgr_post.draws()

array([[[ 1.30511e+02,  9.75941e-01,  7.47178e-01, ...,  2.85270e-03,
          1.01211e-01,  1.14241e-01],
        [ 1.26794e+02,  8.16985e-01,  7.13934e-01, ...,  1.10709e-01,
          3.40024e-02,  8.85645e-02],
        [ 1.30074e+02,  9.94407e-01,  6.84990e-01, ...,  7.31923e-02,
          1.11121e-01,  6.49327e-02],
        [ 1.29947e+02,  9.83667e-01,  7.72327e-01, ...,  3.66501e-02,
          7.89994e-02,  1.18494e-01]],

       [[ 1.30388e+02,  9.89454e-01,  7.47178e-01, ...,  8.69072e-02,
          3.37897e-02,  9.25661e-02],
        [ 1.28357e+02,  1.00000e+00,  7.13934e-01, ...,  1.03536e-01,
          8.78874e-02,  7.51175e-02],
        [ 1.29825e+02,  9.89659e-01,  6.84990e-01, ...,  1.45836e-01,
          1.17892e-01,  9.41959e-02],
        [ 1.30160e+02,  9.85570e-01,  7.72327e-01, ...,  5.90831e-02,
          8.81251e-02,  6.52294e-02]],

       [[ 1.30241e+02,  9.88424e-01,  7.47178e-01, ...,  2.16155e-02,
          7.24129e-02,  1.02635e-01],
        [ 1.28479e+02,  1.00000e+00,  7.13934e-01, ...,  9.62566e-02,
          4.77785e-02,  1.24287e-01],
        [ 1.30376e+02,  1.00000e+00,  6.84990e-01, ...,  9.49858e-02,
          8.78894e-02,  3.76091e-02],
        [ 1.30277e+02,  9.92040e-01,  7.72327e-01, ..., -5.04970e-03,
          8.21757e-02,  7.99845e-02]],

       ...,

       [[ 1.29980e+02,  9.92327e-01,  7.47178e-01, ...,  9.34758e-02,
          8.81039e-02,  8.87766e-02],
        [ 1.27094e+02,  8.98094e-01,  7.13934e-01, ...,  1.67347e-01,
          7.64573e-02,  1.55398e-01],
        [ 1.30398e+02,  9.36709e-01,  6.84990e-01, ...,  7.45936e-02,
          9.61824e-02,  1.78366e-01],
        [ 1.30322e+02,  9.36647e-01,  7.72327e-01, ...,  9.29926e-02,
          1.03869e-01,  1.98546e-02]],

       [[ 1.29567e+02,  9.44852e-01,  7.47178e-01, ...,  7.06839e-02,
          9.29076e-02,  7.62933e-02],
        [ 1.29801e+02,  9.69785e-01,  7.13934e-01, ...,  5.52811e-02,
          1.21212e-01,  2.58968e-02],
        [ 1.29925e+02,  9.67583e-01,  6.84990e-01, ...,  9.35320e-02,
          7.62850e-02,  2.78032e-02],
        [ 1.30498e+02,  1.00000e+00,  7.72327e-01, ...,  1.28179e-01,
          9.27318e-02,  7.85656e-02]],

       [[ 1.29641e+02,  9.26465e-01,  7.47178e-01, ...,  2.57477e-02,
          4.39165e-02,  1.69162e-01],
        [ 1.29784e+02,  8.43336e-01,  7.13934e-01, ...,  6.84627e-02,
          1.58820e-01,  1.42270e-01],
        [ 1.30276e+02,  9.76230e-01,  6.84990e-01, ...,  6.78183e-02,
          7.87947e-02,  8.21316e-02],
        [ 1.30363e+02,  9.83632e-01,  7.72327e-01, ...,  1.22037e-01,
         -2.18041e-03,  8.66952e-02]]])

Code

cgr_post.draws_pd()

         lp__  accept_stat__  ...  cgr_pred[49]  cgr_pred[50]
0     130.511       0.975941  ...      0.101211      0.114241
1     130.388       0.989454  ...      0.033790      0.092566
2     130.241       0.988424  ...      0.072413      0.102635
3     128.701       0.882168  ...      0.075803      0.041847
4     129.255       1.000000  ...      0.128673      0.084683
...       ...            ...  ...           ...           ...
3995  129.794       1.000000  ...      0.037539      0.086167
3996  130.408       1.000000  ...      0.074031      0.071042
3997  130.322       0.936647  ...      0.103869      0.019855
3998  130.498       1.000000  ...      0.092732      0.078566
3999  130.363       0.983632  ...     -0.002180      0.086695

[4000 rows x 59 columns]

Code

cgr_post |> DataFrame

1000×16 DataFrame
  Row │ iteration  chain  μ          σ          lp       n_steps  is_accept  a ⋯
      │ Int64      Int64  Float64    Float64    Float64  Float64  Float64    F ⋯
──────┼─────────────────────────────────────────────────────────────────────────
    1 │       501      1  0.0885972  0.036928   84.5361      3.0        1.0    ⋯
    2 │       502      1  0.0940632  0.0379282  82.8881      1.0        1.0
    3 │       503      1  0.0901949  0.0350495  83.3674      3.0        1.0
    4 │       504      1  0.0847086  0.0315686  82.4833      3.0        1.0
    5 │       505      1  0.0782442  0.0368907  85.4057      3.0        1.0    ⋯
    6 │       506      1  0.081535   0.0384244  85.8435      3.0        1.0
    7 │       507      1  0.0781813  0.0500135  85.2838      3.0        1.0
    8 │       508      1  0.0802804  0.0412565  86.1106      3.0        1.0
  ⋮   │     ⋮        ⋮        ⋮          ⋮         ⋮        ⋮         ⋮        ⋱
  994 │      1494      1  0.0859129  0.0408179  85.6162      7.0        1.0    ⋯
  995 │      1495      1  0.0833438  0.0437737  85.8445      3.0        1.0
  996 │      1496      1  0.0812184  0.0364029  85.3886      3.0        1.0
  997 │      1497      1  0.0812184  0.0364029  85.3886      1.0        1.0
  998 │      1498      1  0.0824591  0.0375846  85.6518      1.0        1.0    ⋯
  999 │      1499      1  0.0648397  0.0572539  83.7585      7.0        1.0
 1000 │      1500      1  0.0717195  0.0492112  85.284       3.0        1.0
                                                  9 columns and 985 rows omitted

taking a look

up and running with PPLs

success

next up: how we can extend this to a robust and helpful workflow

break?

☕

comparisons in `Turing.jl`

Code

n_draws = 1_000; n_chains = 4

# a no U-turn sampler, with 2000 adaptive steps and a target acceptance rate of 0.65
NUTS_sampler = NUTS(2_000, 0.65)

# a Hamiltonian Monte Carlo sampler, with a step size of 0.05 and 10 leapfrog steps
HMC_sampler = HMC(0.05, 10)

# a Metropolis-Hastings sampler, using the default proposal distribution (priors)
MH_sampler = MH()

# a 'compositional' Gibbs sampler (Metropolis within Gibbs) - sampling μ with MH and σ with NUTS
Gibbs_sampler = Gibbs(MH(:μ), NUTS(2_000, 0.65, :σ))

run_mcmc = function(sampler)
    return prepare_data().Δc |> 
      data -> corrosion_growth(data) |>
      model -> sample(model, sampler, MCMCThreads(), n_draws, n_chains)
end

probabilistic programming

overview

from BUGS/JAGS to Stan

from BUGS/JAGS to Stan

from BUGS/JAGS to Stan

PPLs

corrosion rates

corrosion rates

loading data

a corrosion growth rate model

running the model

taking a look

taking a look

up and running with PPLs

break?

comparisons in Turing.jl

from `BUGS`/`JAGS` to `Stan`

from `BUGS`/`JAGS` to `Stan`

from `BUGS`/`JAGS` to `Stan`

comparisons in `Turing.jl`